Chlamee nitrate analysis

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(cache = TRUE)

library(tidyverse)
library(cowplot)
library(broom)
library(readxl)
library(janitor)
library(plotrix)
library(here)
library(growthTools)
library(rootSolve)

Read in data

treatments <- read_excel(here("data-general", "ChlamEE_Treatments_JB.xlsx")) %>%
    clean_names() %>%
    mutate(treatment = ifelse(is.na(treatment), "none", treatment)) %>% 
    filter(population != "cc1629")
nitrate <- read_csv(here("data-processed", "nitrate-abundances-processed.csv"))

this is the step that gets us the growth rate estimates

growth_rates_n_AICc <- nitrate %>%
    filter(population != "COMBO") %>% 
    mutate(ln.fluor = log(RFU)) %>% 
    group_by(well_plate) %>% 
    do(grs = get.growth.rate(x = .$days, y = .$ln.fluor,id = .$well_plate, plot.best.Q = F)) 

Pull out the things we want

growth_sum_n_AICc <- growth_rates_n_AICc %>%
    summarise(well_plate, mu = grs$best.slope, n_obs = grs$best.model.slope.n,
              slope_r2 = grs$best.model.slope.r2,
              best.model_r2 = grs$best.model.rsqr,
              best.model = grs$best.model, best.se = grs$best.se,
              contents = grs$best.model.contents) 

key <- nitrate %>% 
    select(well_plate, population, nitrate_concentration) %>% 
    distinct(well_plate, .keep_all = TRUE)

AICc_growth_rates <- growth_sum_n_AICc %>% 
    select(-contents) %>% 
    mutate(IC_method = "AICc")

AICc_growth_rates2 <- left_join(AICc_growth_rates, key, by = "well_plate")
# write_csv(AICc_growth_rates2, here("data-processed", "nitrate_exp_growth_w_growthtools_AICc.csv"))

exp_params <- growth_sum_n_AICc %>% 
    unnest(contents %>% map(tidy, .id = "number"))  ## pull out the slopes and intercepts etc.

exp_params_aug <- growth_sum_n_AICc %>% 
    unnest(contents %>% map(augment, .id = "number")) %>% ### now add the fitted values
    rename(days = x)

exp_wide <- exp_params %>% 
    spread(key = term, value = estimate)

all_preds <- left_join(exp_params_aug, exp_wide)
all_preds2 <- left_join(all_preds, key, by = "well_plate") %>% 
    group_by(well_plate) %>% 
    mutate(B1 = mean(B1, na.rm = TRUE)) %>% 
    mutate(B2 = mean(B2, na.rm = TRUE)) %>% ### here we do some wrangling to get the colour coding right for our plots
    mutate(exponential = case_when(best.model == "gr" ~ "yes",
                                   best.model == "gr.sat" & days < B1 ~ "yes", 
                                   best.model == "gr.lag" & days < B1 ~ "yes", 
                                   best.model == "gr.lagsat" & days < B2 & days > B1 ~ "yes",
                                   TRUE ~ "no"))

all_preds2 %>%
    ggplot(aes(x = days, y = .fitted, group = well_plate, color = best.model)) + geom_line() +
    geom_point(aes(x = days, y = y, shape = exponential)) +
    facet_grid(nitrate_concentration ~ population) + ylab("Ln(RFU)") +xlab("Days")

Fit the Monod model to the growth rates

monod_fits <- AICc_growth_rates2 %>% 
    mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>% 
    rename(estimate = mu) %>% 
    group_by(population) %>% 
    do(tidy(nls(estimate ~ umax* (nitrate_concentration / (ks+ nitrate_concentration)),
                data= .,  start=list(ks = 1, umax = 1), algorithm="port", lower=list(c=0.01, d=0),
                control = nls.control(maxiter=500, minFactor=1/204800000))))

get the fitted values

prediction_function <- function(df) {
    
    monodcurve<-function(x){
        growth_rate<- (df$umax[[1]] * (x / (df$ks[[1]] +x)))
        growth_rate}
    
    pred <- function(x) {
        y <- monodcurve(x)
    }
    
    x <- seq(0, 1000, by = 0.1)
    
    preds <- sapply(x, pred)
    preds <- data.frame(x, preds) %>% 
        rename(nitrate_concentration.x = x, 
               growth_rate = preds)
}

bs_split <- monod_fits %>% 
    select(population, term, estimate) %>% 
    dplyr::ungroup() %>% 
    spread(key = term, value = estimate) %>%
    split(.$population)


all_preds_n <- bs_split %>% ### here we just use the fitted parameters from the Monod to get the predicted values 
    map_df(prediction_function, .id = "population")


all_predsn_2 <- left_join(all_preds_n, treatments, by = c("population")) %>% 
    distinct(ancestor_id, treatment, nitrate_concentration.x, .keep_all = TRUE)
all_growth_n2 <- left_join(AICc_growth_rates2, treatments)

all_growth_n2 %>% 
    mutate(estimate = mu) %>%
    mutate(nitrate_concentration = as.numeric(nitrate_concentration)) %>% 
    ggplot(aes(x= nitrate_concentration, y= estimate)) + 
    geom_point() +
    geom_line(data=all_predsn_2, aes(x=nitrate_concentration.x, y=growth_rate, color = treatment), size = 1) +
    facet_grid(treatment ~ ancestor_id) +
    geom_hline(yintercept = 0.1, linetype = "dotted") +
    ylab("Exponential growth rate (/day)") + xlab("Nitrate concentration (uM)") 

monod_wide <-  monod_fits %>% 
    select(population, term, estimate) %>% 
    spread(key = term, value = estimate)

m <- 0.1 ## set mortality rate, which we use in the rstar_solve
monod_curve_mortality <- function(nitrate_concentration, umax, ks){
    res <- (umax* (nitrate_concentration / (ks+ nitrate_concentration))) - 0.1
    res
}   

Find R*

rstars <- monod_wide %>% 
    mutate(rstar = uniroot.all(function(x) monod_curve_mortality(x, umax, ks), c(0.0, 50))) %>% ## numerical
    mutate(rstar_solve = ks*m/(umax-m)) ## analytical

rstars2 <- left_join(rstars, treatments, by = "population") %>%
    distinct(population, ks, umax, .keep_all = TRUE)

rstars2 %>% 
    group_by(treatment) %>% 
    summarise_each(funs(mean, std.error), rstar_solve) %>% 
    ggplot(aes(x = reorder(treatment, mean), y = mean)) + geom_point() +
    geom_errorbar(aes(ymin = mean - std.error, ymax = mean + std.error),width = 0.1) +
    ylab("R* (umol N)") + xlab("Selection treatment") + geom_point(aes(x = reorder(treatment, rstar), y = rstar, color = ancestor_id), size = 2, data = rstars2, alpha = 0.5) +
    scale_color_discrete(name = "Ancestor")